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ABSTRACT 

We have developed a new model for analysing light curves of planetary transits when there are 
starspots on the stellar disc. Because the parameter space contains a profusion of local minima 
we developed a new optimisation algorithm which combines the global minimisation power 
of a genetic algorithm and the Bayesian statistical analysis of the Markov chain. With these 
tools we modelled three transit light curves of WASP-19. Two light curves were obtained on 
consecutive nights and contain anomalies which we confirm as being due to the same spot. 
Using these data we measure the star's rotation period and velocity to be 11.76 ± 0.09 d and 
3.88 ± 0.15kms _1 , respectively, at a latitude of 65°. We find that the sky-projected angle 
between the stellar spin axis and the planetary orbital axis is A = 1.0° ± 1.2°, indicating 
axial alignment. Our results are consistent with and more precise than published spectroscopic 
measurements of the Rossiter-McLaughlin effect. 

Key words: planetary systems — stars: fundamental parameters — stars: spots — stars: 
individual: WASP-19 



1 INTRODUCTION 



It was who first put forward the idea of using planetary 

transits to detect starspots, through an anomalous brightening when 
the planet passes over the starspot d uring the transit of the host star. 
As transit surveys such as WASP dPollacco etai]|2006l) and HAT 
(Bakos et al. 2004) detect ever more transiting planets, the number 
of known planets orbiting an active star will incre ase. iRabus et"al] 
d2009h : IPont et alj ( 120071) and IWinn et all feOlObl) have all shown 
how a planet crossing a starspot during transit can create a small 
increase in the received flux from the star. This occurs because the 
spot is generally cooler than the surrounding photosphere, so less 
light is lost when a planet is occulting the spot than when the planet 
is in front of the unspotted parts of the photosphere. One of the 
remaining problems is to model the effects of both the transit and 

spot accurately and precis ely. 

ISanchis-Oieda et all feoill) . iNutzman et all J201lh and 
iDesert etaf (2011) used photometric observations to show that 
it is possible to calculate the obliquity of the system when there 
are light curves of multiple transits affected by the same spot(s). 
A large spi n-orbit misalignment has be en found in this way for 
HAT-P-11 JSanchis-Oieda & Winnll201ll) . For a multiple planetary 
system, starspots have also been used to test the alignment of 
the stellar spin axis agai nst the orbital planes of the planets 
(Sanct uVOieda et al]|2012t) . The obliquity of a planetary system 
helps understand which mechanism was predominant in the 
dynamical evolution of the system iSanchis-Oieda et alj [20121 : 
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IWinnetal.ll2010ah . A low obliquity would follow from the idea 
that the planet formed at a large distance from its host star and, 
through tidal interactions with the protoplanetary disc, suffered 
orbital decay. Larger obliquities are expected when orbital decay 
occurre d due to gravitatio nal interactions from other bodies in the 
system dSchlaufmanll2010l) . 

At present there are three main ways to measure the rotation 
period of a star. The first is to us e photom etric rotational modu- 
lation over many months or years iHail 19721) . The second method 
uses radial velocity mesurements to find the projected rotational ve- 
locity, v sin I, which gives a lower limit on v and thus the rotation 
period. The third method, presented by ISilva- Valid J2008h . is the 
idea of measuring the rotation period of a star by using a transiting 
planet crossing a starspot. This opens up the possibility to allow the 
rotation period of a star to be found from two sets of photometry 
from a 2m-class ground-based telescope. 

Once the rotation period of a star is known it is possible 
to find t he age of the star according to the gyrochronology rela- 
tionship (Barnel [20071 e.g .). It is also possible, when combined 
with v sin/, to find I, the inclination of the stellar spin axis to- 
wards the observer. To find the spin-orbit alignment of a sys- 
tem, iFabrvckv & Wind te009h showed that we require three pa- 
rameters, orbital inclination i, the inclination of the stellar spin 
axis towards the observer I, and the sky-projected spin-orbit 
alignment A. The quantity A can be found by using two dif- 
fere nt methods. The first method is th e Rossiter-McLaughlin ef- 
fect teossitei|[l924 iMcLaughlrnll 19241) . whil e the second method 
uses planetary transits crossing over starspots dSanchis-Oieda et all 
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201 ll; ISanchis-Oieda& Wiiir] l201ll; ISanchis-Oieda et alj 120121 : 



Nutzman et alj201ll ; |Pesert et alfcOllF 

After observing three light curves of WASP- 19 with the aim 
of obtaining accurate physical properties, we discovered that two 
of our datasets contai ned a starspot ano maly, Starspots can affect 
the shape of a transit dSilva-Valioll2010h and if not correctly mod- 
elled can lead to biased measurements of the system parameters. 
To achieve our original goal to obtain precise measurements of the 
system properties we decided to develop a new model capable of 
modelling both the transit and starspots simultaneously. With a pre- 
cise known position of the spot at two close but distinct times we 
would then be able to calculate the obliquity of the system and com- 
pare this to the val ues found from me a surement of the Ro ssiter- 
McLaughlin effect dHellier et alj|201 ll ; lAlbrecht et alj 2012h and 



check if WASP- 19 follows the theory put forward by Winn et al. 
(2010a) that cool stars will have low obliquity (see also iTriaua 
1201 lb . This would also allow us to measure the rotation period of 
the s tar and compare i t to the value found by photometric modula- 
tion drlebb et aljiolol) . 

Section|2] describes the IDlQ code PRISM and how it models 
both the planetary transit and the star-spot. Section[3]describes the 
optimisation algorithm used to find the optimal solution together 
with their associated uncertainties. Section|4] gives an overview of 
the observations and the manner in which they were collected. Sec- 
tion|5]shows the best fitting results of the model for both the system 
and starspots. Section[6] reviews previous work on WASP-19 and 
compares the results of this work. 




Figure 1. An output model of a transit coupled with a starspot using PRISM. 
The transit chord is represented by the two horizontal black lines. The black 
disc to the left is the planet. A dark starspot has been placed on the curved 
stellar surface to show how PRISM projects an elliptical shape onto the stel- 
lar disc for a circular spot. 



2 MODELLING TRANSITS AND STARSPOTS: 
INTRODUCING PRISM 

When dealing with starspot anomalies in transit data a common 
course of action is to model the trans it first and then deal with 
the starspot based on th e residuals (e.g. lSanchis-Oieda et alj[201ll ; 
iMacieiewski et al]|201lh . This can unfortunately lead to unknown 
uncertainties and biases in the measurements of the stellar and 
planetary radii, incl ination and limb darkening coefficients, (e.g. 
iBallerini et all2012h . This is because when a starspot is on the vis- 
ible part of the stellar disc it reduces the received flux by an amount 
AF S pot. When the planet transits the star it blocks A.Fpi anc t of the 
stellar flux. The depth of the transit is the fractional amount of flux 
blocked by the planet, which depends on the ratio of the areas of 
the planet and star. Without a starspot and in the absence of limb 
darkening (LD) the equation for the ratio of the radii is: 



'R P \ 2 _ AF pl 

anct 



F 



(1) 



where F is the total flux of the unspotted star and R p and Rs are the 
radii of the planet and star. If a starspot is placed on the stellar disc 
and causes a decrease in stellar flux of AFspot, the above equation 
becomes: 



AF 



planet 



(F - AF SI 



(2) 



where a is the ratio of the transit depths in the spotted and unspot- 
ted cases. Because AFpot > for a cool spot, the transit gets 



1 The acronym IDL stands for Interactive Data Language and is a trade- 
mark of Exelis Visual Information Solutions. For further details see 

http : //www . ittvis . com/ Product Services/ IDL . aspx 



deeper (a > 1). Neglecting this would result in an incorrect mea- 
surement of the ratio of the radii, 4 E - 

To obtain accurate measurements of the system and spot pa- 
rameters we created an IDL computer code to model both the plan- 
etary transit and starspots on the stellar surface. PRISfvQ (Plane- 
tary Retrospective Integrated Star-spot Model) uses a pixellation 
approach to create the modelled star on a two-dimensional ar- 
ray in Cartesian coordinates (see Fig.Q~|l. This makes it possible 
to model the transit , limb darke ning and starspots on the stellar 
disc simultaneously. ISilvald2003l) used a similar model to describe 
the starspots on HD 209458, but with the drawback of having to 
use fixed system parameters from previous results. PRISM is set to 
use the standard quadratic limb darkening law and uses the frac- 
tional stellar and planetary radii (the radii scaled by the semimajor 
axis, r s , p = R S:P /a). PRISM requires ten parameters to model the 
system:- 

• the ratio of the radii, — = 

• the sum of the fractional radii, r„ + r B — Rp+Rs 



• the linear limb darkening coefficient, ui 

• the quadratic limb darkening coefficient, 112 

• the orbital inclination, i 

• a reference transit midpoint, To 

• the longitude of the centre of the spot, 9. Longitude is defined 
to be zero degrees at the centre of the stellar disc. 

• the latitude of the centre of the spot, tj>. Latitude is defined to 
be zero degrees at the north pole and 180 degrees at the south pole. 

• the spot size, r spo t, in degrees. 

• the spot contrast, p spo t, which is the surface brightness of the 
spot versus the immaculate photosphere 



Available from http://www.astro.keele.ac.uk/~jtr 
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Figure 2. Five light curves showing how the shape of the spot anomaly 
changes with the size of the spot relative to that of the planet. The spot con- 
trast was also modified for each light curve to maintain an approximately 
constant spot anomaly amplitude. This amplitude gives a lower limit for the 
size of the spot. The spot sizes are labelled on the left of the plot and the 
spot constrasts on the right. The 8.5° spot (green solid line) is the repre- 
sentation of when the spot is of equal size to the transiting planet. There 
is a degeneracy between the spot radius and contrast, which can be broken 
when modelling data of sufficiently high precision and cadence. 

To ensure sufficient numerical resolution, the diameter of the planet 
is hard-coded to be 100 pixels, and the size of the star in pixels is 
scaled according to the specified ratio of the radii. When modelling 
a starspot, PRISM projects a circular spot onto the curved surface 
of a star. Because of this it is able to account for spots which are 
visible on the edge of the star even if their centre is beyond the 
limb. 



2.1 Sample light curves 

When a spot anomaly is viewed during a transit the total flux re- 
ceived increases for a dark spot. The total change in fl ux is based 
on the surface area and contrast of the spot ( Silva 2003). Therefore 
there is a degeneracy between these two parameters. Fig.|2] shows 
example light curves for anomalies of approximately the same am- 
plitude and due to a range of spot sizes and constrast. It is possible 
to discern three regimes from this diagram. Firstly, when the spot 
is of a similar size to the planet the shape of the spot-occultation 



is an inverted 'V. This is due to the fact that the amount of time 
the planet spends fully eclipsing the spot is very small compared 
to the duration of the partial eclipse phases. Secondly, for a larger 
spot, both the peak and base of the spot-transit increase, because the 
planet reaches the spot earlier and spends more time fully eclipsing 
the spot. Thirdly, for a smaller spot, the peak broadens due to the 
planet fully eclipsing the spot for longer while the base shortens 
due to the fact that the total duration is shorter. These three distinct 
shapes allow the degeneracy between the spot radius and contrast 
to be broken for data of sufficient precision and time sampling. 

It is also apparent that the amplitude of the spot-transit gives 
a lower limit on the size of the spot, below which the spot is too 
small to give such an amplitude even if its contrast is zero. In Fig.|2] 
the 2.5° spot has a contrast of zero and is still unable to achieve the 
same change in flux as the other spots. 



3 OPTIMISATION ALGORTHMS: INTRODUCING 
GEMC 

Our first attempt at fitting real data with PRISM utilised a Monte 
Carlo Markov Chain (MCMC) algorithm. This was introduced in 
order to use Bayesian methods to find the best fit and associated 
errorbars. We found that the problem with this approach was that 
the many local minima in the parameter space tended to trap our 
MCMC chains, resulting in poor mixing and convergence. This 
could be solved by using a large number of iterations, but such 
an approach was ill-suited to PRISM due to the significant amount 
of processing time required per iteratiorjf]. We found that MCMC 
chains required up to 10 6 iterations to converge properly, depend- 
ing on how often they got stuck in local minima, which equated to 
about a week of calculation time. 

Our solution to this problem was to implement a genetic algo- 
rithm (GA). A GA mimics biological processes by spawning suc- 
cessive generations of solutions based on breeding and mutation 
operators from the previous generation. By performing these oper- 
ations the new solutions are generated based on the fitness of the 
parent solutions, not a perturbation of their parameters. Because 
of this a GA can be considered as a global optimiser where solu- 
tions can jump large distances across the solution space. The effi- 
ciency of a GA at fi nding the global solut ion is d e monst rated by 
ICharbonneau] (l995) but, as discussed by iRajpaull j2012h . it does 
have some limitations. These are primarily that GAs are ill-suited 
to Bayesian statistics, and that they are good at finding where the 
global solution is but poor at locating its exact position. 

Our initial answer to the latter problem was to allow the GA to 
find and constrain the global solution and then to use the MCMC al- 
gorithm to perform the error analysis for this solution. This allowed 
us to reduce the computation time from seven to five days. Dissat- 
isfied with the fact that two different optimisation algorithms had 
to be used, one to locate the global solution and the other to obtain 
parameter uncertainties, we decided to develop a new optimisation 
algorithm, which combined the global optimisation power of the 
GA but also able to perform Baysian statistics on the solutions. We 
call this new algorithm GEMC0 (Genetic Evolution Markov Chain). 
GEMC is based on a Differential E volution Markov Chain (DE-MC) 
put forward bv lTer Braakl ( 120061) . 

3 A single evaluation of a model appropriate for WASP-19, with 70 data- 
points, takes PRISM typically 0.7 s using a 2.7 GHz duel core desktop com- 
puter 

4 Available from http://www.astro.keele.ac.uk/^jtr 
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GEMC begins by randomly generating parameters for N 
chains, within the user-defined parameter space, and then simul- 
taneously evolves the chains for X generations. At each generation 
the chains are evaluated for their fitnesfl The parameters of the 
fittest member undergo a ±1% perturbation and its fitness is then 
re-evaluated. If the fitness has improved it is accepted but if the fit- 
ness has deteriorated it is accepted based on a Gaussian probability 
distribution: 



exp 




(3) 



where (n — 1) is the previous generational chain and n is the current 
generational chain being evaluated. The next step is to then evolve 
the other chains. This is accomplished in a similar way as a GA, in 
that the chain parameters are modified by incorporating the param- 
eters of another chain. But unlike a GA where a member is picked 
by a weighted random number and then the digits of each parame- 
ter are crossed over with the digits from a different member, GEMC 
directly perturbs the parameters of each chain in a vector towards 
the best-fitting chain. The size of this perturbation is between zero 
and twice the distance to the best-fitting chain, allowing the chain 
to not only move towards but also to overshoot the position of the 
best-fitting chain. 

An example would be a two-dimensional function f(x,y). 
The difference between a given chain and the best-fitting chain in 
this case is (Ax, Ay). This difference is then multiplied by a ran- 
dom scalar 7 for each parameter, where 7 is in the interval [0,2], 
and then added to the given chain's parameters (xq, yo) to form the 
new potential solution f(x\,y\). 



xi = x + j x Ax 
yi=yo + lyAy 



(4) 
(5) 



When 7 = the parameter is not perturbed, 7 = 1 the parameter 
equals the current best fitting value and when 7 = 2 the parameter 
is perturbed to the opposite position. This allows the potential solu- 
tion to travel large distances across the parameter space unimpeded 
by local peaks. After the parameters have been perturbed the chain 
is then re-evaluated and is selected using the same method as the 
best fitting chain. 

GEMC runs in two stages. The first stage, called the 'burn in', 
is used to find the optimal solution to the data using the above 
method. After this the second stage starts in which each chain un- 
dertakes an independent MCMC run. The starting points for each 
MCMC chain lie close but not exactly at the optimal solution. In 
essence what we have is the same outcome from running a GA to 
find the best fitting solution and to use this to tightly constrain the 
starting parameter range of an MCMC run. 

When GEMC was used in conjunction with PRISM to find the 
best fitting solution to the same dataset as above, the computational 
time reduced dramatically, from five days to 14 hours using a large 
parameter range (see Section 5). When the parameter range was set 
to the same as used by the GA or the MCMC, GEMC was able 
to produce the best fitting solution and similar uncertainties in the 
fitted parameters as the MCMC within 10 hourfl 



5 A solution's fitness was found by calculating the 1/x 2 value. 
e GEMC is able to produce similar parameter uncertainties as an MCMC 
run with only 1000 iterations (taking 11 min to calculate), but for statisti- 
cal certainty the MCMC section of GEMC was allowed to run for 50 000 
iterations (9.3 hr). 




10 20 30 40 50 
Generation 



Figure 3. The evolution of the fittest chain (solid line) and the mean fitness 
chain (dashed line) from each generation. The maximum peak was found in 
five generations. The fitness is measured as 1 — f(x, y). 



To demonstrate the power of GEMC at finding the optimal so- 
lution of a rugged parameter space we chose to test it against the 
functio n used to test the genetic algorithm PIKAIA bv lCharbonneaj 
d 19951) : 



f(x, y) = [16:r(l — x)y(l — y) sin(n7ra;) sm(nny)] 2 
x,y 6 [0,1], n= 1,2,... 



(6) 



The optimal so lution to this fun ction lies at the centre 
(/(0.5,0.5) = D. ICharbonneaul jl995l) showed that it took PIKAIA 
with a population of 100 solutions up to 20 generations to find the 
global maximum peak but even after 100 generations it still had 
not found the global maximum point, confirming the GA inability 
to find best solutions with precision. While looking at Fig.|4]we can 
clearly see that GEMC, using a population of only 40 solutions, has 
found the global maximum peak within 10 generations and then 
went on to find the global maximum point within 20 generations. 
We can also see from Fig.[3]fhe power of GEMC. The global max- 
imum peak was actually found at the fifth generation and all solu- 
tions were very close to the global maximum point by the twentieth 
generation. This performance indicates that the required burn-in for 
GEMC is extremely short and as such greatly reduces the computing 
time required to find the global solution. 



4 OBSERVATIONS AND DATA REDUCTION 

Three transits of WASP- 19 were observed in February 2010 using 
the 3.6 m New Technology Telescope (NTT) operated at ESO La 
Silla, Chile. The instrument used was EFOSC2, operated in imag- 
ing mode and with a Gunn r filter (ESO filter #786). In this setup 
the CCD covers a field of view of (4.1') 2 with a plate scale of 
0.12" px _1 . No binning or windowing was used, resulting in a dead 
time between consecutive images of 83 s. The exposure time dura- 
tion were 60-90 s. The moon was bright and relatively close to the 
target star. The pointing of the telescope was adjusted to allow five 
good comparison stars to be observed simultaneously with WASP- 
19 itself. We were able to keep the telescope autoguiding through 
all observations. An observing log is given in TableQ] 
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Figure 4. For a comparison with lCharbonnead Jl995l) . OEMC results for A 1 " = 40 chains and for X = 100 generations. The global maximum peak and global 
maximum point have been discovered by the 10th and 20th generations respectively. By the 40th generation all 40 chains have found the global maximum 
point. 
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Table 1. Log of the observations presented in this work. N ohs is the number of observations. 'Moon ilium.' and 'Moon dist.' are the fractional illumination of 
the Moon, and its distance from WASP-19 in degrees, at the midpoint of the transit. 



Date 


Start time 


End time 


N ohs 


Exposure 


Filter 


Airmass 


Moon 


Moon 


Aperture 


Scatter 




(UT) 


(UT) 




time (s) 






ilium. 


dist. 


sizes (px) 


(mmag) 


2010/02/24 


06:18 


09:34 


68 


90 


Gunn r 


1.14 — y 2.30 


0.742 


85.5 


42, 60, 100 


0.573 


2010/02/25 


00:44 


04:26 


76 


60-90 


Gunn r 


1.40 -4 1.04 


0.818 


78.1 


52, 70, 90 


0.464 


2010/02/28 


04:01 


07:41 


74 


90 


Gunn r 


1.04 ->■ 1.42 


0.996 


53.0 


44, 64, 88 


0.499 



These observations were experimental for two reasons. Firstly, 
the NTT is an alt-az telescope fitted with an image derotator. This 
means that the path of light from each star through the telescope 
is continually changing, raising the possibility of correlated noise 
due to any optical imperfections. Secondly, the NTT is fitted with 
an actively controlled thin primary mirror designed to provide the 
best possible focus for normal observing strategies. Defocussing 
such a telescope might lead to a point spread function (PSF) which 
is variable in time, and thus correlated noise via flat-fielding errors. 

In practise we found that, whilst careful attention had to be 
paid to the amount of defocussing, the NTT is perfectly capable of 
producing high-quality light curves whilst a long way out of focus 
due to stable symmetric PSFs. Our observations used this approach 
and ar e not plagued by co rrelated noise. This situation is similar to 
that of Wi nn et all d2009l) . who successfully ob served WASP-4 us- 
ing the Magellan Baade telescope. By contrast, iGillon et id mi) 
encountered serious problems in obtaining photometry of WASP- 
4 and WASP-5 with the ESO Very Large Telescope. This problem 
was attributed to the need to turn off the active optics system in or- 
der to achieve strong defocussing, and our results support the con- 
tention that this is not a general problem with alt-az telescopes or 
active-optics systems. 

We reduce d the data in an identical fashion to 
ISouthworth et alj j2009j|bHcl |2010|) . In short, we performed 
aperture photo metry using an IDL implementation of DAOPHOT 
dStetson|[l987l) . and adjusted the aperture sizes to obtain the best 
results (see TableQJ. A first order polynomial was then fitted to the 
outside-transit data whilst simultaneously optimising the weights 
of the comparison stars. The resulting data have scatters ranging 
from 0.464 to 0.573 mmag per point versus a transit fit using 
PRISM. The timescale used is HJD/UTC. 



5 DATA ANALYSIS 

We began by selecting a search space for each parameter. As dis- 
cussed in Section[3] the ability of GEMC to find the global minima 
in a short amount of computing time meant we were able to search 
a large area of parameter space to avoid the possibility of missing 
the best solution. The parameter search range used in analysing the 
WASP-19 datasets are given in Table|2] 

First we modelled the three datasets of WASP-19 separately 
using PRISM, finding that the modelled parameters were within l-cr 
of each other (Table[2j. We then modelled all three datasets simulta- 
neously. The ensuing parameters agreed with the individual results 
found previously, but we were unable to get as good a fit to the data. 
The reason for this seems to be the LD coefficients, which are in 
comparatively poor agreement when the three light curves are fit- 
ted individually. The scatter around the weighted mean is x2 = 2.2 
for the linear coefficient and 1.9 for the quadratic coefficient. This 
situation co uld be caused by the in fluence of the starspot on the LD 
coefficients. iBallerini et alj d2012l) found that starspots can affect 



LD coefficients by up to 30% in the ultraviolet, with a weaker ef- 
fect expected at redder wavelengths. If we assume a 10% variation 
in the LD coefficients for our r-band data, the coefficients move 
into l-cr agreement between the datasets. 

5.1 Photometric results 

As the combined fit to the three datasets has significantly larger 
residuals than individual fits, we based our final results on the in- 
dividual fits to the data. The final photometric parameters for the 
WASP-19 system are given in Table[3]and are weighted means plus 
l-cr uncertainties of the results from the three individual fits. Fig.[5] 
compares the light curves to the best-fitting models, including the 
residuals. 

The results from modelling the spot anomalies on the nights of 
2010/02/24 and 2010/02/25 confirm that they are due to the same 
spot rotating around the surface of the star, as the spot sizes and 
contrasts are in good agreement. Fig.[6]is a representation of the 
stellar disc, the spot and the transit chord for the two nights of ob- 
servations. 

From the positions of the starspot at the time of the transits 
on the nights of 2010/02/24 and 2010/02/25, it is possible to cal- 
culate the rotational period of the star and the sky-projected spin 
orbit alignment of the system using simple geometry. The spot 
has travelled 24.52° ± 0.28° in 1.015 ± 0.001 orbital periods, 
giving a rotational period of P Iot = 11.76 ± 0.09 d at a lati- 
tude of 65°. Combining this with the stellar radius found below, 
we calculate the latitudinal rotational velocity of the star to be 
w (65°) = 3.88 ± 0.15kms _1 . The positions of the spot finally 
yield a sky-projected spin orbit alignment of A = 1.0° ± 1.2° for 
WASP-19. 

We have collected the a vailable times of mid-transit for 
WASP-19 from the literature fifebb et alj 1201(1 lAnderson et~al] 
l201lMDragomir et al.l20lTI ; lAlbrecht et al.l2012f) . All timings were 
converted to the HJD/TDB timescale and used to obtain a new or- 
bital ephemeris: 

T = HJD/TDB 2 454 775.33754(18) + 0.78883942(33) x E 

where E represents the cycle count with respect to the reference 
epoch and the bracketed quantities represent the uncertainty in the 
final digit of the preceding number. Fig.|7] and Table|4] show the 
residuals of these times against the ephemeris. We find no evidence 
for transit timing variations in the system. 

5.2 Physical properties of the WASP-19 system 

Now we have measured the photometric properties of WASP-19 
we can proceed to the determination of its physical ch aracteristics. 
We un dertook this analysis following the method of Southworth 
(2009), which uses the parameters measured from the light curves 
and spectra, plus tabulated predictions of theoretical models. We 
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Table 2. Derived photometric parameters from each lightcurve, plus the interval within which the best fit was searched for using GEMC. 



Parameter 


Symbol 


Search interval 




ZUlU/UZ/Zj 


ZUiU/UZ/Zo 


Radius ratio 


r p /r s 


0.05 to 0.30 


0.1435 ± 0.0014 


0.1417 ±0.0013 


0.1430 ±0.0008 


Sum of fractional radii 


r s + r p 


0.10 to 0.50 


0.3298 ± 0.0041 


0.3300 ± 0.0025 


0.3311 ±0.0044 


Linear LD coefficient 


Ml 


0.0 to 1.0 


0.314 ± 0.095 


0.501 ± 0.083 


0.438 ± 0.077 


Quadratic LD coefficient 


U2 


0.0 to 1.0 


0.192 ±0.023 


0.222 ± 0.019 


0.226 ± 0.009 


Inclination (degrees) 


i 


70.0 to 90.0 


78.97 ±0.39 


78.92 ± 0.37 


78.91 ± 0.44 


Transit epoch (HJD/UTC) 




±0.5 in phase 


2455251.79628 ±0.00014 


2455252.58506 ± 0.00010 


2455255.74045 ± 0.00012 


Longitude of spot (degrees) 


e 


-90 to +90 


-9.54 ±0.15 


14.98 ±0.13 




Latitude of Spot (degrees) 


4> 


0.0 to 90.0 


64.93 ± 0.32 


65.37 ± 0.21 




Spot angular radius (degrees) 


f'spot 


0.0 to 30.0 


15.01 ±0.21 


15.18 ±0.15 




Spot contrast 


pspot 


0.0 to 1.0 


0.777 ± 0.011 


0.760 ± 0.017 





0.0005 ■ 



tj 0.0000 



w -0.0005 — 



-0.00101 



200 



400 600 
Cycle number 



800 



1000 



Figure 7. Residuals of the available times of mid-transit versus the orbital ephemeris found in this work. The three timings from this work are the cluster of 
three points around cycle number 600. 



Table 3. Combined system and spot parameters. 



Parameter 


Symbol 


Value 


Radius ratio 


r p /r s 


0.1428 ± 0.0006 


Sum of fractional radii 


r s + r p 


0.3301 ± 0.0019 


Linear LD coefficient 


Ml 


0.427 ± 0.049 


Quadratic LD coefficient 


U2 


0.222 ± 0.008 


Inclination (degrees) 


i 


78.94 ± 0.23 


Spot angular radius (degrees) 


^spot 


15.13 ±0.12 


Spot contrast 


Pspot 


0.771 ± 0.010 


Stellar rotation period (d) 


Prot 


11.76 ±0.09 


Projected spin orbit alignment (degrees) 


A 


1.0 ± 1.2 



adopted the values of i, r p /r 3 and r s + r v from Table[5] and the 
stellar properties of effective temperature T c ff = 5440 ± 60 K 
(Maxted e tai1l201ll) . velocity amplitude K s = 257 ± 3ms 
(Helli eretal. | l201lh and metal abundance [§] = 0.02 ± 0.09 
(Helli eretal.ll201ll) . 

An initial value of the velocity amplitude of the planet, K p , 
was used to calculate the physical properties of the s ystem using 
standa rd formulae and the physical constants listed bv lSouthworthl 
fcoill) . The mass and [^] of the star were then used to obtain the 
expected T c g and radius, by interpolation within a set of tabulated 
predictions from stellar theoretical models. K p was iteratively re- 
fined until the best agreement was found between the observed and 
expected T c ff , and the measured r B and expected — . This was per- 
formed for ages ranging from the zero-age to the terminal-age main 
sequence, in steps of 0.01 Gyr. The overall best fit was found, yield- 
ing estimates of the system parameters and the evolutionary age of 
the star. 

This procedure was performed se parately using five different 
sets of stellar theoretical models (see Southworth1 l2010l) . and the 



Table 4. Times of minimum light of WASP-19 and their residuals versus 

the ephemeris deri ved in this w o rk. 

References: (1) iHebbetalJ J2010h ; (2) lAlbrecht et al.1 <2012l) ; (3) 
lAnderson etld] 1201 ll) ; (4) This work; (5) lDragomir et alj fepl ll)~ 



Time of minimum 


Cycle 


Residual 


Reference 


(HJD/TDB - 2400000) 


no. 


(HJD) 




54775.33757 ±0.00020 


0.0 


0.00004 


1 


55168.96839 ±0.00011 


499.0 


-0.00001 


2 


55183.16748 ±0.00007 


517.0 


-0.00003 


3 


55251.79657 ±0.00014 


604.0 


0.00003 


4 


55252.58544 ±0.00010 


605.0 


0.00005 


4 


55255.74077 ±0.00012 


609.0 


0.00003 


4 


55580.74238 ± 0.00058 


1021.0 


-0.00020 


5 



Table 5. Physical properties of the WASP-19 system. 



Parameter 


Value 




Stellar mass ( Mq) 


0.904 ± 0.040 ±0 


021 


Stellar radius ( R©) 


1.004 ± 0.016 ±0 


008 


Stellar surface gravity (cgs) 


4.391 ± 0.008 ±0 


003 


Stellar density ( pq ) 


0.893 ±0.015 




Planet mass ( Mj up ) 


1.1 14 ± 0.036 ±0 


017 


Planet radius ( Rj up ) 


1.395 ± 0.023 ±0 


011 


Planet surface gravity (ms~ 2 ) 


14.19 ±0.26 




Planet density (pj up ) 


0.384 ± 0.011 ±0 


003 


Equilibrium temperature 


2067 ± 23 




Safronov number 


0.02852 ± 0.00057 ± 


00023 


Semimajor axis (AU) 


0.01616 ± 0.00024 ±0 


00013 


Age (Gyr) 


n .+2.7+0.7 
L1 -° -2.3 -1.5 
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Figure 5. Transit light curves and the best-fitting models. The residuls are 
displayed at the base of the figure. 




Figure 6. Representation of the stellar disc, starspot and transit chord for 
the two datasets containing spot anomalies. 

spread of values for each output parameter was used to assign a sys- 
tematic error. Statistical errors were propagated using a perturba- 
tion algorithm. An alternative set of physical properties was calcu- 
lated using a calibration of ste llar properties base d on well-studied 
eclipsing binary star systems (Enoch et al. 2010), with calibration 
coefficients from South worth] d201 1 ). 

The final results of this process are in reasonable agreement 
with themselves and with published results for WASP- 19. The fi- 
nal physical properties are given in Table.[5]and incorporate sepa- 
rate statistical and systematic errorbars for those parameters which 



depend on the theoretical models. The final statistical errorbar for 
each parameter is the largest of the individual ones from the solu- 
tions using each of the five different stellar models. The systematic 
errorbar is the largest difference between the mean and the individ- 
ual values of the parameter from the five solutions. One point to 
note is that the inferred age of the star is rather large, particularly 
given its rotation period and activity level. The age is governed pri- 
marily by the input T e g and [^],soa check of these spectroscopic 
parameters would be useful. 



6 SUMMARY AND DISCUSSION 

We have introduced the PRISM code to model a planetary transit 
over a spotted star, and the optimisation algorithm GEMC for find- 
ing the global best fit and associated uncertainties. While GEMC 
still requires significant computing time to calculate parameter un- 
certainties via Markov chains, the speed at which it can find the 
optimal solution is a large improvement over the long burn-in cur- 
rently required by an MCMC routine. 

We have applied PRISM and GEMC to three transit light curves 
of the WASP- 19 planetary system. Two of the light curves are of 
consecutive transits and show anomalies due to the occultation of 
a starspot by the planet. The measured latitudes and longitudes of 
the spot during the two transits were used to calculate the rotation 
period of the star and the sky-projected obliquity of the system. Our 
model assumes that the spot anomaly can be represented by a circu- 
lar spot of uniform brightness. It is quite likely that the "spot" is in 
fact a group of smaller spots with lower contrasts, but investigation 
of this puts extreme demands on data quality and quantity which 
are practially impossible to satisfy for ground-based observations. 

We find a rotation pe riod of P m t = 11.76 ± 0.09 d at 
a latitude of 65°, whereas iHebb et al.l feoich found a P Iot of 
10.5 ± 0.2 d from rotational modulation of the star's brightness 
over several years. The latter value comes from the spot activity 
over the whole visible surface of the star, whereas our value is 
for a specific latitude. The difference betw een these two numbers 
may therefore indicate differential rotation. I Anderson et al. U20H 
used the measured CaH&K line activity index, logR' HK , to in- 
fer P ro t = 12.3 ± 1.5 d using the activity-rotation calibration by 
Mamajek & Hillenbrand ( 2008), which i s in good agreeme nt with 
the values measured by ourselves and bv lHebb et alj ( 120100 . 

We find a rotational velocity of t>(65°) = 3.88 ± 0.15 km s~ 
for WASP- 19 A, which in the absence of differential rotation 
would yield an equatorial rotatio n velocity of ^(90°) = 4.30 ± 
0.15kms _1 . iHellier et al.1 J20T 1) reported a spectroscopic mea- 
surement for vsin/ of 5.0 ± 0.3kms _1 and assumed this value 
represented the equatorial velocity. They included it as a prior when 
modelling the Rossiter-McLaughlin effect, finding a final value of 
11 sin/ = 4.6 ± 0.3kms _1 . This last measurement is appropriate 
for the latitude at which the planet transits, and may differ from 
ours due to the effect of starspots on radial velocity measurements 
taken during transit. 

We find a sky-projected obliquity of A = 1.0° ± 1.2° for 
WASP- 19, which is in agreement with but more precise than pub- 
lished values based on observations of the Rossiter-M cLaughlin ef- 
fect ( 4.6° ± 5.2°. lHellieret alj|201ll ; 15° ± 11°, lAlbrecht et al.l 
|2012|) . A give s the lower boundary of t he true spin-orbit angle, tp. 
As stated by iFabrvckv & Wiiinl {2009), finding a small value for 
A can be interpreted in different ways. The spot method could al- 
low us to determine ifr, rather than just A, given light curves of 
three or more transits all showing anomalies due to the same spot. 
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But with only two light curves it is difficult to be sure that tj) lies 
close to A. We calculated a minimum rotation period of WASP-19 
of 5.5 d, for the extreme case that the orbital axis was aligned with 
the line of sight. This result disagrees w ith previous measurements 
faebb et al.ll20lol ; I Anderson et al]|201ll) . Whilst we are unable to 
determine the true value of tp with the data in hand, we find no evi- 
dence for a spin-orbit misalignment in the WASP-19 system. With 
a low obliq uity and coo l host s tar, WASP-19 follows the idea put 
forward by IWinn et alj fcOlOah that planetary systems with cool 
stars will have a low obliquity. It also lends weight to the idea that 
WASP- 19 b formed at a much greater distance from host star and 
suffered orbital decay through tidal interactions with the protoplan- 
etary disc. 
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